-
Notifications
You must be signed in to change notification settings - Fork 173
Simplify interpolate
#4582
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Simplify interpolate
#4582
Conversation
92645d1
to
aa30bf8
Compare
expr_args = extract_arguments(ufl.as_ufl(expr)) | ||
is_adjoint = len(expr_args) and expr_args[0].number() == 0 | ||
v = Argument(v.dual(), 1 if is_adjoint else 0) | ||
V = Argument(V.dual(), 1 if is_adjoint else 0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Put 0 if it is not there, otherwise keep increasing the count
V = Argument(V.dual(), 1 if is_adjoint else 0) | |
arg_numbers = [a.number() for a in expr_args] | |
number = 0 if len(expr_args) == 0 or min(arg_numbers) > 0 else max(arg_numbers) + 1 | |
V = Argument(V.dual(), number) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is_adjoint
goes away.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Similar fixes will be needed in UFL
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In UFL we're checking for contiguous argument numbers, I think this approach is more fool-proof: just get the smallest number that's missing
V = Argument(V.dual(), 1 if is_adjoint else 0) | |
arg_numbers = set(a.number() for a in expr_args) | |
number = min(set(range(len(expr_args) + 1)) - arg_numbers) | |
V = Argument(V.dual(), number) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
expr_args
can only contain zero or one arguments so I think this logic is more confusing than is_adjoint
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We should then assert len(expr_args) <= 1
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we don't need to fix this. We have lots more issues if we want to support rank 3 tensors, even symbolically.
class Interpolate(ufl.Interpolate): | ||
|
||
def __init__(self, expr, v, | ||
def __init__(self, expr, V, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is very unrelated, but I think that a much more friendly interface is to allow either or both left and right arguments to be a primal FunctionSpace.
Right now we do this under the hood
Interpolate(Function(V1), V2) -> Interpolate(Function(V1), Argument(V2.dual(), 0))
It'd be reasonable to have a similar shortcut for the adjoint. When the left argument is a FunctionSpace, we would then automatically create the Argument for it.
Interpolate(V1, Cofunction(V2.dual())) -> Interpolate(Argument(V1, 0), Cofunction(V2.dual()))
And supplying two FunctionSpaces is a perfectly natural interface:
Interpolate(V1, V2) -> Interpolate(Argument(V1, 1), Argument(V2.dual(), 0))
Of course we need to arbitrarily decide who gets the lowest number, the more intuitive numbering that produces the forward Interpolation is to go from right to left.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thoughts @dham ?
10fb533
to
f8c2318
Compare
test -> trial fix hypre-ads
cofunction docstring
fix
f8c2318
to
31f1cec
Compare
if isinstance(V, functionspaceimpl.WithGeometry): | ||
# Need to create a Firedrake Argument so that it has a .function_space() method | ||
expr_arg_numbers = {arg.number() for arg in extract_arguments(expr) if not is_dual(arg)} | ||
is_adjoint = len(expr_arg_numbers) and expr_arg_numbers == {0} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is_adjoint = len(expr_arg_numbers) and expr_arg_numbers == {0} | |
is_adjoint = expr_arg_numbers == {0} |
v = Argument(v.dual(), 1 if is_adjoint else 0) | ||
if isinstance(V, functionspaceimpl.WithGeometry): | ||
# Need to create a Firedrake Argument so that it has a .function_space() method | ||
expr_arg_numbers = {arg.number() for arg in extract_arguments(expr) if not is_dual(arg)} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are we sure that dropping all the dual arguments the right thing? We should only drop the dual argument of the source space (if it appears in the expression)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We need to verify/ensure that extract_arguments
returns free arguments only. extract_free_arguments
this would when visiting a FormBase, drop the dual argument in position 0.
v = Argument(v.dual(), 1 if is_adjoint else 0) | ||
if isinstance(V, functionspaceimpl.WithGeometry): | ||
# Need to create a Firedrake Argument so that it has a .function_space() method | ||
expr_arg_numbers = {arg.number() for arg in extract_arguments(expr) if not is_dual(arg)} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We need to verify/ensure that extract_arguments
returns free arguments only. extract_free_arguments
this would when visiting a FormBase, drop the dual argument in position 0.
expr_args = extract_arguments(ufl.as_ufl(expr)) | ||
is_adjoint = len(expr_args) and expr_args[0].number() == 0 | ||
v = Argument(v.dual(), 1 if is_adjoint else 0) | ||
V = Argument(V.dual(), 1 if is_adjoint else 0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we don't need to fix this. We have lots more issues if we want to support rank 3 tensors, even symbolically.
Simplify
interpolate
andInterpolate
. This removes the argument renumbering 0 -> 1 in the primal expression if V is a FunctionSpace. Hence this is an API change and should be merged after the next release. PR #4572 warns users about this.